Relativistic quantum effects of Dirac particles simulated by ultracold atoms 
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Quantum simulation is a powerful tool to study a variety of problems in physics, ranging from 
high-energy physics to condensed-matter physics. In this article, we review the recent theoretical 
and experimental progress in quantum simulation of Dirac equation with tunable parameters by 
using ultracold neutral atoms trapped in optical lattices or subject to light-induced synthetic gauge 
fields. The effective theories for the quasiparticles become relativistic under certain conditions in 
these systems, making them ideal platforms for studying the exotic relativistic effects. We focus on 
the realization of one, two, and three dimensional Dirac equations as well as the detection of some 
relativistic effects, including particularly the well-known Zitterbewegung effect and Klein tunneling. 
The realization of quantum anomalous Hall effects is also briefly discussed. 
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I. INTRODUCTION 

As first pinpointed by Richard Feynman in 1982, a 
quantum computer could be used to simulate a quan- 
tum system efficiently without experiencing an exponen- 
tial explosion and slowdown, which a classic computer 
would presumably encounter when simulating quantum 
phenomena [1]. Although quantum computers are not 
yet available, one can still hope to create analog quan- 
tum simulators which are designed to have the same 
Hamiltonian as another system of interest. The recent 
rapid progress in quantum coherent control [2] of neu- 
tral atoms, photons, or ions in some systems makes them 
become ideal quantum simulators [3] , which include neu- 
tral atoms in optical lattices (OLs), arrays of cavities, 
trapped ions, quantum dots, superconducting circuits, 
nuclear magnetic resonance (NMR) and so on. 

Analog quantum simulators would be able to pro- 
vide highly controllable platforms to study difficult prob- 
lems in various fields, ranging from high-energy physics, 
condensed-matter physics to quantum chemistry. Some 
recent notable achievements in quantum simulation are 
listed as follows: addressing neutral atoms in OLs to 
simulate the Bose-Hubbard model with the quantum 
phase transition from a superfluid to a Mott insula- 
tor [4-6] , as well as to simulate the antiferromagnetic spin 
chains [7]; simulating a quantum magnet [8], frustrated 
Ising model [9], relativistic dynamics [10,11], and open- 
system dynamics [12] with trapped ions; implementing a 
hydrogen molecule quantum simulation with NMR [13]. 
Besides, many other experimental proposals of quantum 
simulation with different simulators are suggested (see 
the review [3] and the references therein). 

One of the most promising quantum simulators is the 
cold atom system since the realization of condensation 



of bosonic atoms [14,15] and fcrmionic atom pairs [16,17] 
has led to num.erous advances with applications in coher- 
ent manipulation of neutral atoms, such as creating cold 
molecules in quantum chemistry (see the reviews [18,19] 
and the references therein). Ultracold neutral atoms in 
OLs are quite suited for mimicking condensed-matter 
physics [20,21], and may shew new light on strongly cor- 
related state problems like high-temperature supercon- 
ductivity [22], fractional quantum Hall effect [23,24], and 
the quantum magnetism [25]. Ultracold atomic gases in 
disordered optical potentials pave the way towards the 
realization of versatile quantum simulators for the in- 
vestigations of Anderson localization [26,27], which may 
help solve some open questions relying on the interplay 
of disorder and interactions [28,29]. In addition, some 
phenomena and effects in high-energy physics and cos- 
mology can be mimicked, such as quantum simulation of 
black holes [30] and cosmic inflation [31] by using Bosc- 
Einstein condensates (BECs), or superstrings [32] and 
supersymmetry [33,34] by using Bose-Fermi mixtures. 

Since the relativistic Dirac fermions were found in 
graphene [35,36], a substantial amount of efforts has 
been devoted to the understanding of exotic relativistic 
effects in solid state systems and the searching for other 
relativistic systems such as the recent discovery of topo- 
logical insulators [37]. Inspired by those exciting results 
and the state-of-the-art technologies in quantum control 
of atoms, one specific topic among ultracold neutral atom 
simulators has arisen, that is, quantum simulation of the 
Dirac equation and the related relativistic effects. 

The Dirac equation ihdt^ = (cq;-p + toc^/3)^, where a 
and /3 are the Dirac matrices, and c is the speed of light, 
was first proposed by Dirac who was seeking a relativistic 
and quantum mechanical equation to describe spin 1/2 
particles with mass m. This equation naturally provides 
a description of the electron spin, which is an assumption 
in the Schrodinger equation. More surprisingly, it pre- 
dicts the existence of antiparticles [38], which was soon 
confirmed by the experimental discovery of positrons. 
The Dirac equation also predicts some exotic effects, and 
the most famous ones are Zitterbewegung (ZB) [39], an 
unexpected trembling motion of free relativistic particles, 
and Klein tunneling (KT) [40] , which describes relativis- 
tic particles penetrating through high and wide potential 
barriers without exponential damping expected in non- 
relativistic tunneling processes. Despite these relativistic 
effects have attracted lots of interest over years, they are 
failed to be directly tested by elementary particles due 
to currently unaccessible experimental techniques. Tak- 
ing a free electron as an example, the frequency and the 
amplitude of ZB is on the order of 10^^ Hz and 10^'^ A, 
respectively, and the electric field gradients requirement 
for observing KT is on the order of 10^^ V/m. All those 
conditions are out of reach in current technologies. 

However, we will see that ultracold atoms loaded in 
some optical lattices or subject to certain synthetic gauge 
fields behave as the relativistic particles under suitable 
conditions. Furthermore, cold atom simulators offer us 



rather more degrees of freedom to control the relativistic 
quasiparticles, such as the dimensions, the effective mass 
and the effective speed of light. Associated with the con- 
trollable atomic interactions and disorder, these systems 
provide us an ideal platform to investigate the interest- 
ing relativistic effects such as the never-before-seen ZB 
and KT for free particles. Although the original Dirac 
equation is specific to fermions with spin 1/2, not bosons 
with an integer spin, we will see that one can still simulate 
the Dirac equation with bosons: a pseudo-spin 1/2 can 
be introduced to the bonons and the dynamics of such 
pseudo-spin should be described by the Dirac equation. 
In this article, we review the recent theoretical inves- 
tigation and experimental proposals in quantum simula- 
tion of the Dirac equation and some related effects by 
using ultracold neutral atoms. This review is organized 
as follows. In Section II, we introduce the realization 
and detection of the two and three-dimensional relativis- 
tic quasiparticles by using cold atoms in OLs with a wide 
range of structures. The realization of quantum anoma- 
lous Hall effect (QAHE) is also briefly introduced. In 
Section III, we introduce another kind of proposals for 
simulating the tunable Dirac equation that are based on 
generating effective gauge fields on bulk atomic gases. 
The simulated Dirac equations do not need OLs and thus 
operates in a continuous regime. In Section IV, we show 
that one can observe the well-known ZB and KT by us- 
ing cold atoms with the schemes reviewed in Section HI. 
Furthermore, we show that an exotic macroscopic KT 
described by the nonlinear Dirac equation (NLDE) ex- 
hibits. Finally, in Section V, we give our summary and 
conclusion of this review, together with some prospects 
on the quantum simulation of relativistic effects and re- 
lated fields in physics. 



II. SIMULATION OF DIRAC EQUATION WITH 
ULTRACOLD ATOMS IN OL SYSTEMS 



In this section, we first review that ultracold neutral 
fcrmionic atoms in a honeycomb (hexagonal) OL, similar 
with electrons in graphene, can mimic 2D massless and 
massive Dirac fermions [41-47]. Then we demonstrate 
that, it is also possible to simulate 2D Dirac fermions in 
a so-called Ta OL [48,49], a line-centered-square (LCS) 
OL [50] , and even in a simple square OL combing with a 
synthetic gauge potential [51-58]. Furthermore, we show 
that 3D Dirac fermions may be realized with cold atoms 
properly loaded in some cubic OLs [59-61]. The real- 
ization of quantum anomalous Hall effect (QAHE) as- 
sociated to the parity anomaly of 2D Dirac fermions is 
also briefiy introduced. In words, the emergent relativis- 
tic Dirac fermions in these OL systems originate from 
the single-particle band dispersion, which provides the 
Dirac cones in the Brillouin zone. Some lattices with 
symmetric structures like hexagonal, Ta and LCS in the 
real space, may directly support the needed structures 
of the momentum space; the simple square lattice may 



fail, and thus needs the aid of external fields. The near- 
half-filling condition allows the Fermi level being close to 
the Dirac points, and in the long wavelength limitation, 
the dynamics of the cold atoms are described well by the 
relativistic theories. The OLs allow one to tune the tun- 
neling probabilities from site to site as well as the atomic 
interactions, and may make the relativistic quasiparticles 
realizable a.nd controllable in these svstems. 
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Fig. 1 The honeycomb OLs. (a),(b) The contours with three 
potentials described in Eq. (1). The minima of the potentials are 
denoted by the solid dots. All V^ are the same in (a), and V^ = 
Vj* = O.QlVg^ in (b). (c) Decomposition of the hexagonal lattice 
as two triangular sublattices A and B with anisotropic tunnelling, 
(d) The Brillouin zone of the hexagonal lattice. The dispersion 
relations are shown in (e) for /3 = 1 (gapless state) and (f) for 
/3 = 2.5 (gapped state). Reproduced from Ref. [41], Copyright 
© 2007 the American Physical Society. 



A. Simulation of 2D relativistic Dirac fermions 

The graphene material, formed with a single layer 
of carbon atoms [35,36] with its emergent massless 
Dirac fermions, has recently attracted strong interest in 
condensed-matter physics. One topic naturally arises is 
to mimic the graphene and the relativistic quasiparticles 
with cold atom in a similar 2D hexagonal lattice [41-47]. 
In addition to the graphene-type lattice, it is also inter- 
esting to search other proposals for quantum simulation 
of 2D Dirac equation with cold atoms in OLs of other 
structures [48-58]. 



1. Mimic graphene: Cold atoms m a honeycomb OL 

Simulating Dirac equations with cold atoms loaded in a 
honeycomb OL was proposed in the first time in Ref. [41]. 



We review this proposal in this section. Considering 
single-component fermionic atoms (e.g., spin-polarized 
atoms ''"K, ^Li, etc.) in a two-dimensional {x-y plane) 
hexagonal OL formed with three detuned standing-wave 
lasers [62] 



Vix,y)= J2 ^"sin^ 

J = l,2,3 



fci {x COS 9j + y sin Oj) + — , (1) 



where 6*1 — 7r/3, 62 = 27r/3, 6*3 = 0, and /c^ is the optical 
wave vector. It is easy to tune the potential barriers Vj* 
by varying the laser intensities along different directions 
to form a standard hexagonal lattice for Vi — V2 — V^, 
and a hexagonal lattice but with a finite anisotropy for 
different Vj* as shown in Fig. 1(a) and 1(b), respectively. 
A hexagonal lattice consists of two sublattices denoted 
by A and i3, as shown in Fig. 1(c). For single-component 
fermionic atoms, the atomic collisions are negligible at 
low temperatures. The tight-binding Hamiltonian of the 
system is then given by 






(aj6,+H.c.), 



(2) 



where {i,j) represents the neighboring sites, a^ and bj 
denote the fermionic mode operators for the sublattices 
A and B, respectively. The tunneling rates Uj depend on 
the tunneling directions in an anisotropic hexagonal lat- 
tice, and we denote them as ti, t2, ^3 corresponding to the 
three different directions as shown in Fig. 1(c). For sim- 
plicity, we assume ti = t2 = t and fs = j3t with /3 being 
the anisotropy parameter. As the atomic tunneling rate 
in an OL is exponentially sensitive to the potential bar- 
rier, this control provides an effective method to control 
the anisotropy of the atomic tunneling by laser intensi- 
ties. The first Brillouin zone of this system has also a 
hexagonal shape in the momentum space with only two 
of the six corners in Fig. 1(d) are inequivalent, corre- 
sponding to two different sites A and B in each cell in 
the real hexagonal lattice, usually denoted as K and K' . 
One can choose K = (27r/a)(l/\/3, 1) and K' = — K, 
where a — 27r/(v3fcL) is the lattice spacing. Taking a 
Fourier transform a\ = (l/yiV) ^j^ exp(ik • Ai)a\^ and 
6] = {\/^/N)Y^^^ex]i{i'k ■ Bj)aj{^, where A,; (Bj) repre- 
sents the position of the site in sublattice A (B) and N 
is the number of sites of the sublattice, the Hamiltonian 
(2) can be diagonalized and the eigenvalues have the ex- 
pression [41] 



J 



Ek = ±tj2 + 13^ + 2 cos{kya) + 4/3 cos{V3k^a/2) cos{kya/2). 



(3) 



As plotted in Fig. 1(e) and 1(f), there are two branches 
of the dispersion relation, corresponding to the ± sign in 
Eq. (3). When < /3 < 2, the two branches touch each 
other, and around the touching points there appears a 
Dirac cone structure. One has the same standard Dirac 
cones as the graphene material with (3 — 1 [35,63,64], 
and the cones squeeze in the x or y direction as /3 de- 
viates from 1, but they still touch each other. When 
/3 > 2, a finite energy gap Ag — \t\{(3 — 2) appears be- 
tween the two branch. So, across the point /? = 2, the 
topology of the Fermi surface changes, corresponding to 
a quantum phase transition without any usual symme- 
try breaking [65]. Such topological phase transition as- 
sociated with pair production (annihilation) events has 
been investigated in Ref. [66] . The evolution of the Dirac 
points in the hexagonal lattice by varying the asymmetry 
hopping and the resulting phase transition was also stud- 
ied in Ref. [44]. With this phase transition, the system 
changes its behavior from a semimetal to an insulator 
at the half filling case (means one atom per cell; note 
that each cell has two sites). Around the half filling, the 
Fermi surface is close to the touching points, and one 
can expand the momentum k around one of the touching 
points K = (fcS, fcO) as (^2, k"y) = (fcO + q,, fcO + Qy). Up 
to the second order of qx and qy, the dispersion relation 
(3) becomes 



E„ 



± 



^l+vlql+vlql, 



(4) 



0, Vx = V3/3ta/2, and Vy == ta^l - /3^/4: 

|i|(/3 - 2), Vx = ^^3/3/2, and 

1 for /3 > 2. This simplified dispersion 



where Ag 
for < /3 < 2; A 
Vy = ta^/Jj2 

relation _Eq is actually a good approximation (named as 
long wavelength approximation) as long as g^, qy ^ l/2a. 
Compared with the standard energy-momentum relation 
for the relativistic Dirac particles, here Ag and v^.y take 
the meaning of rest energy and the velocity of light 
respectively. The wave function for the quasiparticles 
around the half filling then satisfies the Dirac equation 
ifidt^ = "Hd^, where the relativistic Hamiltonian Hd is 
given by 

Hd ^ VxCTxPx +VyayPy + Agffz, (5) 

where 'Jx,y,z are the three Pauli matrices. 

Through an analogy to the graphene physics, we have 
shown that by controlling the lattice anisotropy, one can 
realize both massive and massless Dirac fermions and 
observe the phase transition between them. This pro- 
posal was recently proved to be experimentally feasible 
in Ref. [45] , where the temperature requirement and crit- 



ical imperfections in the laser configuration are consid- 
ered in detail. Even in the presence of a harmonic con- 
fining potential, the Dirac points are also found to sur- 
vive [67]. In the presence of atomic interactions, the 
many-body physics of Dirac particles in graphene-type 
lattices, such as novel BCS-BEC crossover [42], topo- 
logical phase transition between gapless and gapped su- 
perfiuid [47] and even charge and bond ordered states 
with p-orbital band of lattices [43,68], have been inves- 
tigated. Notably, the realization of ultracold quantum 
gases in a hexagonal OL was reported in a very recent 
experiment [69] , which paves the important way to mimic 
the relativistic Dirac fermions and the aforementioned 
beyond-graphene physics with controllable systems. 



2. Pseudospin-1 massless Dirac fermions: Atoms in a 
Ti, /line- centered- square OL 

It is interesting to note that cold atoms in OLs with 
other structures of symmetries can also simulate the 2D 
relativistic Dirac fermions. One of such examples was 
proposed in Ref. [48]. It showed that atoms trapped in 
the 7i3 lattice can behave as the massless Dirac fermions 
with pseudospin S' = 1, instead oi S — 1/2 for those in 
the hexagonal lattice. The so-call Ts lattice, illustrated 
in Fig. 2(a), has a unit cell with three different lattice 
sites, one sixfold coordinated site H and two threefold 
coordinated sites A and B. 



The Ta OL can be experimentally realized through 
three counterpropagating pairs of laser beams with the 
same wavelength 3/2a and linearly polarized with the 
electrical field in the x-y plane, which are similar with the 
laser setup for creating a hexagonal OL. Given a polariza- 
tion of a pair of lasers on the y axis, the other two pairs 
are obtained by rotating 27r/3 around the z axis. The 
Schrodinger equation for the Tjs OL filled with fermionic 
atoms in the tight-binding approximation is given by 

S*H(RH) = -t^*A(RH + T,) + *B(Rff+r,+i),(6a) 

3 

s*„(R«) = -iX!*ff(^"-^j)' ae{A,B}. (6b) 



Here ^q,(Rq) is the amplitude of the wave function on 
sublattice a {— A,H,B), and the tj connect nearest 
neighbors. Solving the Eq. (6), one can obtain the 
energy-momentum relationship 



J 



So(k)=0, 

£;±(k) == ±t^& + 4{cos[(vi - V2) • k] + cos[vi • k] + cos[v2 • k}] 

I 



(7a) 
(7b) 



where E± exhibit two Dirac points, K and K' as seen in Fig. 2(b), and Eq represents a flat band. For 



the near half-filling case, one arranges the component 
a = {A, H, B} into a pseudospin triplet and expands the 
momentum k in the vicinity of K, and thus obtains the 
effective Hamiltonian of Dirac-Weyl form given by 



n 



DW 



vpS ■ p, 



(8) 



which describes the massless Dirac fermions with total 
pseudospin 5 = 1. Here vp = 3ta/^/2 is the Fermi veloc- 
ity, p is the momentum operator in the xy plane, and the 
pseudospin vector operator S = {SxTSy,Sz) iSx,y,z are 
3x3 matrices which satisfy angular momentum commu- 
tation relations) reflected the three inequivalent lattice 
sites per unit cell in Ta- A comparison on the topologi- 
cal properties between the hexagonal lattice and the Ta 
lattice was given in Ref. [49], where the former is shown 
to be topological insulator and the latter to be trivial 
insulator in the framework of quantum spin Hall states. 







Fig. 2 (a) The Ts lattice. It is characterized by translation vec- 
tors VI = (3/2; -V3/2)a and V2 = (3/2; V3/2)a. Open circles 
mark the two sublattices A and B, forming a hexagonal lattice. 
Solid circles mark the hub sites H forming a (larger) triangular 
lattice, (b) Contour plot of E^ (k) in units of the hopping energy 
i, cf. Eq. (7b). The dashed hexagon defines the first Brillouin 
zone, K = 27ra-l(l/3; -v/3/9) and K' = 27ra-l(l/3; ^3/9) are 
two nonequivalent Dirac points. Reproduced from Ref. [48], Copy- 
right © 2009 the American Physical Society. 
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Fig. 3 (a) Schematic illustration of the LCS lattice, energy con- 
tour of the optical potential in Eq. (9) with (b) Vi = 2V2 and (c) 
Vi = 2.2V2, and dispersion relation on the LCS lattice with (d) 
A' = and (e) A' = t. Reproduced from Ref. [50], Copyright 
© 2010 the American Physical Society. 

Another example of quantum simulation of 2D mass- 
less Dirac fermion was proposed to be realized in a LCS 



OL [50]. The massless Dirac fermions there also have 
pseudospin- 1, and may exhibit a perfect all-angle KT (we 
will discuss KT in the Section IV), i.e., the barrier is com- 
pletely transparent for all incident angles [50]. The LCS 
lattice is schematically illustrated in Fig. 3(a), with three 
sublattices denoted by the black (sites A), red (sites B), 
and green (sites C) points, respectively. The LCS lat- 
tice can be realized in the ultracold atomic system by 
applying six detuned standing-wave laser beams, four of 
which are applied along the e^ and e^ directions with op- 
tical wave vectors kL and 2kL, respectively; the other two 
are applied along the (e^, -I- ej,)/v2 directions with opti- 
cal wave vector 2kL and relative phase 7r/2, respectively. 
Thus the whole optical potential is 



V{x,y) = Vi[ain^ (kLx) + sin"^ (kLy) +sm'^{2kLx) 



+ sm^(2kLy)] + V2{sm'^ kL(x + y) + - 



kL{x - y) + -^} 



(9) 



with tunable potential amplitudes Vi and V2, and lat- 
tice constant a — Tr/fc^. Two examples of energy con- 
tours of LCS OL are plotted in Fig. 3(b) and 3(c), in 
which the potential minima are marked with the black 
and blue points. For Vi — 2V2, V{x,y) is the same at 
sites A-C so that their site energies ea = es = ec- While 
for Vi 7^ 2V2, es = ec 7^ ^A and one may set eb = 
— e^ — A' for the symmetry consideration. The Hamil- 



^i C, Ci 



tE 



{^■J)^i^3 



tonian of LCS lattice 'Hlcs = X] 
can be diagonalized with three branches through a sim- 
ilar procedure for Eq. (2). One is a flat band with en- 
ergy Eq 

E±=± 



A' and the other two dispersive bands are 

A^2 + 4i2[cos2(fc^a/2) + cos^{kya/2)]. Typical 

energy bands are shown in Fig. 3(d) for A^ = and Fig. 
3(e) for a;, = t. 

At A' — 0, only one nonequivalent Dirac point appears 
at K = {n/a, ir/a). In the vicinity of the Dirac cone, the 
ultracold atoms behave as the massless Dirac fermions 
with pseudospin-1, described by Eq. (8) with vp = ta 
in this case. Interestingly, the massless Dirac fermions 
in the HSC lattice are topologically different from those 
in the hexagonal and Ta OLs presented previously be- 
cause there are not two nonequivalent Dirac cones but a 
single one in the first Brillouin zone. Note that the solid- 
state material with a single (or odd number) Dirac cone, 
named as strong topological insulator, has attract consid- 
erable attention recently [37]. Besides, the flat band with 
Dirac cones meeting at the same energy, which means the 
particle can have an infinite effective mass (flat band) or a 
zero effective mass (Dirac fermions), may result in inter- 
esting atom dynamics in OLs, such as different behaviors 
of atomic localization for bosons and fermions [70] . 



3. Dirac fermions in a square OL with a gauge field 

The massless Dirac fermions can also be simulated by 
using ultracold atoms in a simple square OL subjected 
to certain light-induced synthetic gauge fields (which will 
be addressed in Section III). The synthetic fields are nec- 
essary because the energy bands of a single square lattice 
do not exhibit the Dirac cone structure. 
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Fig. 4 (a) The atomic levels and the interactions between atoms 
and laser fields, (b) Schematic representation of the experimental 
setup with the two laser beams incident on the cloud of atoms, 
(c) Schematic of the square optical lattice and the designed phase 
factor (denoted by arrows) . (d) The scheme of overlapping the two 
state-selective optical lattices. Reproduced from Ref. [52], Copy- 
right © 2009 the American Physical Society. 

The first scheme along this line was proposed in 
Ref. [52], where the fermionic atoms are loaded in a 
square OL subjected to a U(l) gauge field generated 
by laser-atom interactions. As shown in Fig. 4(a) and 
(b), the noninteracting atoms of mass m with four lev- 
els in the x-z plane are coupled with three laser beams. 
The ground state |1) is coupled to the excited state |3) 
via a laser field with the corresponding Rabi frequency 
riie"*''^^ and the state |2) is coupled to the excited state 
1 3) via a laser field with the corresponding Rabi fre- 
quency ri2e*'^^^. The cold atoms are trapped in two state- 
selective optical potentials as shown in Fig. 4(c) and (d), 
say atoms with internal states |1), |2), and |3) are trapped 
in odd columns (sublattice A), while atoms with inter- 
nal states 1 4) are trapped in even columns (sublattice 
B). The two sublattices make up a 2D rectangular lat- 
tice with the lattice spacings a^ and a^, especially a 2D 
square lattice for ax = a-z- In the basis {|1), |2), |3), |4)}, 
the total Hamiltonian of this system can be written as 



ns = Ho + h 



where Hq 



( Va i7ie^9i^ n^ \ 

Va ^26-"^^^ 

f^ie-*?!^ rJae"?^^ Va 

\ fl3 Vb / 



, (10) 



are the two state-selective periodic potentials. Diagonal- 
izing the Hamiltonian (10) yields two degenerate dark 
states (with eigenenergy £■ = 0) 



Ixi) =cos0|l)-sin6'e*«^|2), 
1X2) = |4), 



(11a) 
(lib) 



and two bright states, where q = qi + q2 and tan 6 = 
|r2i|/|ri2|. If the motions of the atoms that initially pre- 
pared in the dark state subspace satisfy the adiabatic 
condition, one can safely use the adiabatic approxima- 
tion and reduce the Hamiltonian to 



2^^t 



H = / d-'r^l 



-^(-zW-A)2+l/^ 
2m 



$1 



+ I d'^r^l 



^2 



$4 



+hne 



V(41 



$1 



$|$4 



(12) 



in the second quantized form, using {|xi),|X2)} as the 
basis. Here $i(r) and $|(r) are field operators corre- 
sponding to annihilating and creating an atom with the 
internal quantum state \i) at coordinate position r re- 
spectively, He = ri3Cos6' and the U(l) adiabatic gauge 
potential A = fi(7 sin ^e^. Taking the tight-binding limit, 
one can superpose the Bloch states to get Wannier func- 
tions Wa(r — r;) and W{,(r — rj) for sublattice A and B, 
respectively. Expanding the field operator in the lowest 
band Wannier functions as 

$i(r)= Y^ am,„e^-''o""'^-''''wa(r-r™„), (13a) 

7n{odd) ,71 

$4(r) = ^ bjn,nWb{r - r™„), (13b) 

m{even).n 

and then substituting them into Eq. (12), one can rewrite 
the Hamiltonian as 

tl = —t 2_^ "m+l,n+l"'"+l + ^ ^m,n+l^rn,n 



[ra{odd) ji] 



"2am rfim+l.n + H.C., 



(14) 



-V^ is kinetic energy operator, and Va,b 



where 7 = 2'Khaz sin is the phase resulted from the 
adiabatic gauge potential. Here the isotropic atomic 
tunneling is assumed, which can be realized by well 
adjusting the intensities of laser beams. Diagonalizing 
the Hamiltonian (14) with 7 = tt and ax = a-z = o,, 
one obtains the quasiparticle energy spectrum E{h) = 
zt2ty^cos'^{kxa) + cos^(fcza), which exhibits two inequiv- 
alent Dirac points in the first Brillouin zone with one 
at K = {n /2a, IT /2a). In the vicinity of the Dirac point 
K, the atoms behave like the massless Dirac fermions 
described by the effective Dirac-like Hamiltonian 

'H = hvF{pxa-x+PzO'z), (15) 

where vp = 2ta/h is the Fermi velocity of this system. 



The 2D massless Dirac fermions also emerge when ul- 
tracold fermionic atoms are trapped in a square OL sub- 
jected to a U(2) non-Abehan gauge field [56], that is to 
say, the phase factor 7 in Eq. (14) is replaced by a 2 x 2 
matrix. This kind of non-Abelian OLs could be created 
by laser-assisted tunneling for atoms in optical superlat- 
tices [71-73]. They are characterized by a constant Wil- 
son loops, and the single-particle spectrum may depict 
a complex structure termed as the Hofstadter "butter- 
fly" [74]. In the proposal, the half- integer quantum Hall 
effect that has been observed in the graphene [35,36] and 
the squeezed Landau vacuum due to anisotropic of exter- 
nal gauge filed are also investigated. 

The external gauge filed can even be replaced by a 
time-dependent optical potential [53,54], which acts to 
shake the lattice, leading to the modification of effec- 
tive tunneling strength. Such driven tunneling can be 
used to generate an artificial staggered magnetic field for 
atoms in a two-dimensional square optical lattice. For 
bosonic gas in this optical system, the zero-temperature 
phase diagram was shown to exhibit a finite-momentum 
superfluid phase with a quantized staggered rotational 
fiux [53]. For fermionic gas, 2D massless Dirac fermions 
may also be provided [54]; and mixing both bosonic and 
fermionic atoms, one may even reach the strongly in- 
teracting regime for 2D Dirac fermions [55], which al- 
lows to investigate the transition between the Dirac and 
the Fermi liquid, and some correspondences with high- 
temperature cuprates. 

Another model with a different time-dependent opti- 
cal potential, which is also used to create an artificial 
magnetic field in cold atoms loaded on a square OL, was 
recently presented in Ref. [58]. The effective dynamics of 
low energy quasiparticles in that system can be described 
by a Dirac Hamiltonian for massless fermions with the 
unusual property that chiral symmetry is broken in the 
tunneling energy term, leading to splitting the doubly de- 
generate Dirac cones into two with tunable slopes. The 
two slopes describe two speeds of light for the relativistic 
quasiparticles. 




Fig. 5 A crystal cell of the edge-centered cubic lattice and the 
axis system adapted in this paper. One crystal cell includes four 
inequivalent lattice sites labeled by 0, 1, 2 and 3, respectively. Re- 
produced from Ref. [59], Copyright © 2010 the American Physical 
Society. 



B. Simulation of 3D relativistic Dirac fermions 

The above mentioned work are all limited in 2D Dirac 
fermions, and the 3D ones in the OLs would be valuable 
to explore since the 3D Dirac fermions differ from 2D ones 
by transport properties, localization, etc. We review in 
the following that cold atoms loaded in an edge-centered 
cubic (ECC) OL can be used to simulate the 3D Dirac- 
like fermions [59]. 

The ECC lattice as shown in Fig. 5 is the three di- 
mensional counterpart of the LCS. In certain conditions 
(established by Eq. (19) below), the cold atoms on the 
lattice sites can be regarded approximately as the ground 
states of anisotropic 3D harmonic oscillators with energy 

4'^ = 5^E.=.,,,.^^ where u;^^ = [9^y/m]p is the 
frequency of the oscillator on site i along the axis v with 
potential l/(see Ref. [59] for its expression). Substituting 
V into the expression of w* , yields 



e(") - E 



-Vai — 4a2 + 4a3 



.(1) 



Er 



\Joix + 



4a3 + -\/-ai 



4a2 + 4a3 



(3) ^ (2) ^ (1) 
S ^9 9 ' 



(16a) 



(16b) 
(16c) 



where ai = Ai/Er{i — 1,2,3) are dimensionless pa- 
rameters with Er = h? {2tt / a)^ / 2m being the recoil en- 
ergy. The 3D harmonic potentials on sites and 1 
are lifted by the bottom energies 1^(0,0,0) = 6A2 and 
V{a/2, 0, 0) = Ai + 2^2 on site and site 1, respectively. 
The ground state energies of these sites are 






(17a) 
(17b) 



Since the ground state energies on lattice sites 1, 2 and 
3 are the same, we may address only one of them, such 
as Eg ' in Eq. (17b). We limit ourself to the energy 
band formed by the ground states of the lattice sites. 
To this end, one has to tune the parameters Ai,A2 and 
A3 (or ai,a2 and as) carefully to satisfy the following 
conditions 



I40) 
1.(0) . 



-£;«!« 2 min[ef,eW], 
41)1 « 2^,. 



(18a) 
(18b) 



The parameter Vj is the well depth that can be estimated 
as Vj = [2y(0, 0, a/4) - V{0, 0, 0) - V{0, 0, a/2)]/2 = A3. 
Equation 19(a) ensures that other energy levels of the os- 
cillators on each sites are separated sufficiently far away 
and then the energy bands are formed by the ground 
states, and Eq. 18(b) guarantees that the calculated 
ground levels on all sites are much lower than the well 
depth so as to local levels can be formed and the ground 
energies can be evaluated by Eq. (16a)-(17b). One can 
estimate that the parameter region with V1/V2 ~ 4 and 
Vs J^ 30-E,- fulfills the above condition. 
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Under the above harmonic approximate conditions, the 
tight-binding Hamihonian of the system in momentum 
space can be worked out as 



Hk = 2t 



cosF^ 

cosFj, 

V cos r^ 



cos Vx COS Vy cos Fj 
(5 

(5 

(5 



(19) 



in the basis {|0), |1), |2), |3)} with |i) being the single par- 
ticle ground state on site i. Here the mid-point of the 
energy on sites and 1 is chosen as the energy reference, 

X, y, z) along 



^x,y,z = ^x,y,zO^I'^ with wave vector ki (z 



- f7(i) 



(0) 



the direction i^2b — E\' ~ E\' vs, the energy difference 
between site 1 and site 0, and t is the hopping constant. 
Diagonalizing -ff/j, one obtains the four band dispersions 
given by 

Ei=E2= (5, (20a) 

E± = ±v/4i2(cos2 F^ -t- cos2 F^ -t- cos^ F^) -I- ,52. (20b) 

The first two dispersions are flat bands stick to the bot- 
tom of the upper band or the top of the lower band, 
depending on the sign of 5. Note that 3D bulk flat band 
has not attracted much attention, in contract to previ- 
ous work based on ID or 2D model [50,56,70]. The last 
two dispersions possess the characteristic of Dirac parti- 
cles near the only one Dirac point K = (tt/a, vr/a, tt/o) at 
the corner of the first Brillouin zone, and they touch each 
other when 5 ^{). At the vicinity of the Dirac point, the 
dispersion (21b) is approximately rewritten as 



E± =±^r 



- p^c*^, 



(21) 



where p — hk is the momentum with k ~ {k'^ + k'^ + klY^'^ 
being the amplitude of 3D wave-vector, c* — at/h is the 
effective light speed, and m* = h^5/{atY is the effective 
rest mass. The dispersion describes 3D massive Dirac 
fermions, and reduced to massless Dirac dispersion, i.e., 
E± — ±c*p, when (5 = 0. 

Apart from this scheme, it has been shown [60] that the 
3D relativistic fermions can also be simulated with ultra- 
cold fermionic atoms in a 3D optical superlattice as an 
extension system of 2D non- Abelian OL in [56] . Interest- 
ingly, the 3D relativistic fermions simulated in this way 
are the so-call naive Dirac fermions in lattice gauge theo- 
ries [75,76], and thus provide a realization of the fcrmion 
doubling problem: Two massless Dirac fermions appear 
in the Brillouin zone, each of them has a different chi- 
rality, corresponding to the right and left-handed parti- 
cles. Furthermore, by tuning the laser intensities, the 
system may allow to decouple the doublers from a sin- 
gle Dirac fermion through inverting their effective mass, 
i.e., a quantum simulation of Wilson fermions [77]. In 
this regime, the atomic gas corresponds to a phase of 
matter, 3D topological insulators, where Maxwell elec- 
trodynamics is replaced by axion electrodynamics [37,78] . 
The chirality symmetry breaking via mass terms is also 
crucial for simulating the Wess-Zumino supersymmetry 



model [34]. 

The 3D Dirac fermions using ultracold atoms are also 
proposed to be realized in a cubic optical lattice sub- 
jected to a synthetic frustrating magnetic field B — 
7r(^o(l,l,l) with (po = h/2ma^ [61]. The tight-binding 
Hamiltonian of the system can be described by 



(id) 



de 



(22) 



where the phase factor is Aij = ^ Ji' A- dl, with the vec- 
tor A — 7r(0, x— J/, y—x) since B — rotA. Diagonalization 
of Hamiltonian (22) yields the dispersion as 



E± = ±2t 



i/(cos2 



(23) 



which exhibit two inequivalent Dirac points again. 
Around the Dirac points, the atoms behave like the 3D 
massless Dirac fermions. It is also shown to be possible 
to give mass to the Dirac fermions by coupling the ultra- 
cold atoms to a Bragg pulse, and to obtain a dimensional 
crossover from 3D to 2D Dirac fermions by varying the 
anisotropy of the lattice. 





0.08 




l^t 






0.06 


\ 


\ ^'^ / 


0.04 




\ / 


0.02 

-8 




\,/ 


2 


-0.1 0.1 0. 



trap edge 



trap center 



|i/t 



Fig. 6 The number <iensity of atoms n per unit cell of the hexago- 
nal lattice as a function of the chemical potential fi (corresponding 
to a re-scaled atomic density profile in a trap) for (a) /9 = 1, and 
(b) (3 = 2.5. A plateau with a width 2/3 — 4 appears for the latter 
case which corresponds to the case when the chemical potential 
sweeps inside the energy gap. (c) The derivative dn/d/i as a func- 
tion of the chemical potential fi for /3 = 1. (d) An enlarged part of 
dn/dfj, at the vicinity of /.t = 0. The linearity of the curve shows the 
linear dispersion relation for the quasiparticles. Reproduced from 
Ref. [41], Copyright © 2007 the American Physical Society. 
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Fig. 7 The number of atoms n per unit cell in Ts OLs at zero 
temperature as a function of fi/t. Inset: density of states, dn/dfi, 
versus fi/t. Reproduced from Ref. [48], Copyright © 2009 the 
American Physical Society. 



C. Detection of the Dirac quasiparticles 

We have shown that cold atoms trapped in OLs with 
a wide range of structures can be used to simulate the 
Dirac equation, and the quasiparticles in these systems 
behave as the relativistic ones. Another important is- 
sue, however, is to experimentally verify the existence of 
these relativistic quasiparticles. The widely used detec- 
tion technique based on the transport measurements for 
the condensed-matter materials is typically not available 
for the atoms, and nevertheless, there are some specific 
detection methods for the trapped atomic gas. For sim- 
plicity, we focus on the aforementioned hexagonal OL 
system [41], and in the following we review two different 
methods to confirm the relativistic quasiparticles and the 
phase transition between massless and massive ones [41]: 
the density profile measurement [79] and the Bragg spec- 
troscopy [80]. However, the detection methods can be 
extended to other systems without loss of generality. 

Density proGle measurement. — The density profile of 
the trapped atoms can be measured through the time- 
of-fiight imaging with the light absorption [79]. Free 
fermions expand with ballistic motion and from the final 
measured absorption images, one can reconstruct the ini- 
tial real-space density profile of the trapped gas. Under 
the local density approximation, the local chemical po- 
tential varies with the radial coordinate by /i = /ip — y(r), 
where (Xq is the chemical potential at the trap center and 
V{r) — ma;^r^/2 is the global harmonic trapping poten- 
tial. At temperature T, the atomic density is given by 



n(p) 



1 

Sf) 



J K'^xj ^yj fJ-jdkxdk. 



y 



(24) 



where 5*0 — STT^/y/ia^ is the area of the first Bril- 
louin zone of the hexagonal lattice, and f{kx,ky,ji) — 
l/{cxp[{Ei,-n)/T] + l} {Ei, inEq.(3)) is the Fcrmi-Dirac 
distribution. At low temperature (T ^ 0), this density 
profile n{fi) is shown in Fig. 6(a) and 6(b) for the pa- 
rameters with massless and massive Dirac quasiparticles, 
respectively. Clearly, a plateau at the atom density n = 1 
in the density profile appears for the gapped phase with 
massive Dirac fermions; but there is no such a plateau for 
the case with massless Dirac fermions. So the plateau is 
associated with massive quasiparticles, and its emergence 
provides an unambiguous signal for the quantum phase 
transition between the two cases. Furthermore, the lin- 
ear dispersion relation for the massless Dirac fermions 
can be confirmed by the derivative dn/dfi since one has 
§^ = :j5§- \Sfi\ around the Dirac cone ioi (3 — 1 and 
Vx — Vy — V. As shown in Fig. 6(c) and (d), ^ is 
linearly proportional to 6fi indeed at the vicinity of the 
half filling. So experimentally, from the measured den- 
sity profile n(y^), one can determine its slope, which sig- 
nals the linear dispersion relation, to confirm the massless 
Dirac fermions. 

This method is also available for the detection of mass- 
less Dirac fermion in Ts OLs [48] . The similar results are 



shown in Fig. 7, but with a sharp jump in the atomic 
density at fx — due to the contribution from the highly 
degenerate flat band. 

Bragg spectroscopy. — The Bragg spectroscopy can pro- 
vide an alternative and complementary method to con- 
firm the linear dispersion relation for the massless Dirac 
fermions and the energy gap for the massive ones [41]. 
In Bragg spectroscopy [80], one shines two laser beams 
on the atomic gas as shown in Fig. 8(a). By fixing the 
angle between the two beams (thus fixing the relative mo- 
mentum transfer q = kj — ki , where k^ denotes the wave 
vector of each laser beam) , one can measure the atomic 
(or photonic) transition rate by scanning the laser fre- 
quency difference uj = uj2 — i^i- From the Fermi's golden 
rule, this transition rate basically measures the following 
dynamical structure factor [80] 

^(q,c^)= Y, |(/kji^s|*ki)|'<5[fi^-i?/k,+i?.kj,(25) 

ki,k2 

where Hb = Eki.ka ^e*i'"|iki)(/k2l + h.c. is the light- 
atom interaction Hamiltonian, and |ikj^) and l/ks) denote 
the initial and the final atomic states with the energies 
Siki and -E/k2 ^^id the momenta ki and k2, respectively. 
At the half filling, the excitations are dominantly around 
the touching point, and we can use the approximate dis- 
persion relation in Eq. (4). For the isotropic case (/? = 1) 
with massless Dirac Fermions, S{q, lu) has the expression 



S{q,uj) 



0, 



{uj < UJr) 
(CJ > LOr) 



(26) 



where LUr — qvp/h (q = |q|) and qr = huj/vp. This 
dynamical structure factor is shown in Fig. 8(b). Note 
that in this case, the lower cutoff frequency ujr is lin- 
early proportional to the momentum difference q, and 
UJr vanishes when q tends to zero. The ratio between 
UJr and q gives the Fermi velocity vp, an important pa- 
rameter as the analogy of the light velocity for conven- 
tional relativistic particles. For the anisotropic case with 
/? > 2, the spectrum in Eq. (4) becomes quadratic with 
E Ki ±(A -I- h^q^/2m,x -\- h^qy/2m,y) for small momentum 
transfer q, where the effective mass m^^.y — f^lS.jv^ . 
The dynamical structure factor in this non-relativistic 
limit becomes 



S{q,uj) = 



0, 



(lu < ui^'y) 



(27) 



where ujp'^ = 2A -I- h^q^ y/'irnx.y- Its form is shown in 
Fig. 8(b). The lower cutoff frequency w^'^ in this case 
does not vanish as the momentum transfer goes to zero. 
This distinctive difference between the dynamical struc- 
ture factors in Eqs. (26) and (27) can be used to distin- 
guish the cases with massive or massless Dirac fermions. 
From the variation of the cutoff frequency uj^'^ as a func- 
tion of the momentum transfer qx,y, one can also experi- 
mentally figure out the important parameters such as the 
energy gap A and the effective masses ttIx and my. 
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Fig. 8 Schematic of the Bragg scattering, (a) Atoms are illumi- 
nated by two laser beams with wave vectors ki , k2 and frequencies 
ui, u)2, respectively, (b) The dynamic structure factors S{q,ui) 
(not scaled) for the massless Dirac fermions (solid line) and for the 
massive ones in the non-relativistic limit (dotted line). The experi- 
mentally measurable quantities iJr and uic'^ give important param- 
eters for the quasiparticles. Reproduced from Ref. [41], Copyright 
© 2007 the American Physical Society. 



D. Quantum anomalous Hall effects 

As a unique electronic behavior of single layer graphene, 
the relativistic Dirac fermions carrying charge in this 2D 
material may result in novel transport phenomena. One 
of the most exotic examples is that, when subjected to 
a perpendicular magnetic field, the so-call half-integer 
quantum Hall effect (QHE) as a result of the uncon- 
ventional relativistic Landau level can be observed in 
graphene [35,36]. For cold atomic systems, such half- 
integer QHE is principally able to be simulated since the 
relativistic particles emerge [56,61], and the additional 
requirement is an effective uniform magnetic field, which 
can be generated by rotation or optical dressing (will be 
discussed in Section HI. A). 

In contrast to the conventional QHE and the half- 
integer QHE, there is another kind of QHE, called as 
quantum anomalous Hall effect (QAHE) because the net 
magnetic flux is zero in each unit cell and there are no 
Landau levels. However, all the three need the breaking 
of time-reversal symmetry. Twenty years ago, Haldane 
proposed a toy model to illustrate such QAHE as a re- 
sult of the parity anomaly of the 2D Dirac fermions in a 
hexagonal lattice [81]. Nevertheless, it is extremely hard 
to realize the Haldane's model experimentally in ordi- 
nary condensed matter systems including the graphene 
because of the unusual staggered magnetic flux assumed 
in the model. Other schemes have been proposed to re- 
alize QAHE in semiconductor systems with topological 
band structures [82-84], where the time- reversal symme- 
try is broken by the spin polarization. Besides, it is possi- 
ble to realize QAHE in graphene in the presence of both 
Rashba spin-orbit coupling and an exchange fleld [85], 
or in a nanopatterned 2D electron gas [86], or in a thin 
film of a three-dimensional topological insulator with an 
exchange field [87]. The QAHE states are promising for 
future device applications because of its potential real- 
ization of dissipationless charge transport. 



Interestingly, several proposals have been suggested to 
realize and detect the QAHE by using ultracold atoms. 
An ingenious scheme is proposed to realize the Haldane' 
model by using cold atoms trapped in a hexagonal OL 
with an artificial staggered magnetic field (Berry curva- 
ture) with hexagonal symmetry in Ref. [88] . The wanted 
external field can be achieved by coupling atomic internal 
state with the standing waves of laser beams. While it 
remains a challenge to experimentally realize such a stag- 
gered magnetic field. In Ref. [89] , it was showed that the 
QAHE can be simulated with the p-orbital band in the 
hexagonal OL through orbital angular momentum polar- 
ization by applying a technique of rotating each lattice 
site around its own center [90]. An expanded version of 
this proposal can be found in Ref. [91]. Besides, it was 
recently demonstrated to realize the QAHE with cold 
atoms trapped in an anisotropic square OL which can be 
generated from available experimental setups of double- 
well lattices with minor modifications [92] . The detection 
of the QAHE in cold atomic systems can be achieved by 
direct detection of the Chern number with the typical 
density-profile- measurement technique [88], or alterna- 
tively determining the topological phase transition from 
usual insulating phase to quantum anomalous Hall phase 
by light Bragg scattering of the edge and bulk states [92] . 
All of these proposals involve experimental techniques of 
implementing cold atoms to be developed. 



III. SIMULATION OF DIRAC EQUATION 

WITH ULTRACOLD ATOMS IN SYNTHETIC 

NON-ABELIAN GAUGE FILEDS 



In the models discussed in Sec. II, the emergency of rel- 
ativistic Dirac fermions (equation) in OL systems gener- 
ally requires the long wavelength approximation, which 
supports the continuous limitation. In this section, we 
introduce another kind of proposals to simulate the tun- 
able Dirac equation that are based on generating effective 
gauge fields [93-101] on bulk atomic gases, which do not 
need OLs and thus operates in a continuous regime [102- 
105]. We will see that the effective relativistic dynamics 
emerge again in the cold atom systems by using atom- 
light interaction to generate effective gauge potentials 
acting on these neutral atoms. Besides, there is no filling 
requirement in these continuous systems, and thus the 
relativistic Dirac dynamics are equivalent for both non- 
interacting bosonic and fermionic in the single-particle 
physics. Moreover, it even holds for the condensation 
of bosons, which provides us a macroscopic counterpart 
of relativistic particles described by the nonlinear Dirac 
equation (NLDE). We first introduce the schemes and 
experiment progress on creating synthetic gauge poten- 
tials in neutral atoms. Then we use two specific systems 
to show that a liner and nonlinear Dirac equation with 
tunable parameters through laser fields can be simulated. 
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A. Synthetic gauge fields for neutral atoms 



Schrodinger equation 



We first briefly review the generation of syntlretic gauge 
fields for neutral cold atoms loaded in OLs. In Section 
II. A. 3, we have mentioned that a U(l) Abelian gauge po- 
tential can be generated by atom-laser interaction [52] or 
temporal modulation of the lattice potential [53,54] for 
atoms trapped in a square OL. This external gauge po- 
tential behaves as a vertical staggered or uniform (when 
7 is spatially independent) magnetic field acting on a 
charged particle moving in a square potential, and thus 
contributes a scalar Peierls phase factor to the atomic 
hopping process. A more straightforward way to gener- 
ate an Abelian gauge potential in OLs can be achieved 
by using laser-assisted tunneling for atoms trapped in 
state-dependent lattice [71,73], where atoms also acquire 
a phase when tunneling between the two neighbor sub- 
lattices. The phase factor can extend to a 2 x 2 matrix, 
corresponding to a square lattice subjected to a U(2) non- 
Abelian gauge potential [94] . The basic idea to create a 
U(N) non- Abelian gauge potential in an OL is to addi- 
tionally consider atomic species with N internal quasi- 
degenerate sub- levels in both two sublattices [72], and 
thus instead of a simple phase, the laser-assisted tun- 
neling would induced a rotation in the internal space, 
described by a iV x A^ matrix F with non-commuting 
component Fx,Fy,Fz. 

An alternative method to create synthetic gauge field 
is rotating the neutral atomic gas trapped in harmonic 
potentials or OLs [106]. This scheme is limited to rota- 
tionally symmetric configurations, and furthermore, can 
not extend to non- Abelian case without additional laser 
fields [107]. Other schemes that can be used to create a 
gauge field include using light beams with orbital angular 
momentum [93-95] or spatially shift laser beams [96-98] 
to interact with atomic gases. Two reviews on the syn- 
thetic gauge potentials for neutral atoms are given in 
Refs. [99] and [108]. 

The relativistic Dirac Hamiltonian always includes a 
spin-orbit coupling term (as shown in Eq. (5)), which 
can be induced from a U(2) non- Abelian gauge potential. 
To achieve this goal without OLs, we can consider the 
scheme based on the adiabatic motion of multiple-level 
atoms in laser fields that generate (near-)degenerate dark 
states [100]. This idea goes back to the work by Wilczek 
and Zee [101], who have shown that non- Abelian gauge 
fields can arise in the adiabatic evolution of quantum sys- 
tems with multiple degenerate eigenstates. We consider 
the adiabatic motion of atoms with N + 1 internal lev- 
els in stationary laser fields, as shown in Fig. 9(a). For 
fixed position r the internal Hamiltonian including the 
laser-atom interaction can be diagonalized to yield a set 
of iV -I- 1 dressed states |Xn(r)) with eigenvalues e„(r), 
where n = l,2,...,A^ + l. The full quantum state of the 
atom describing both internal and motional degrees of 
freedom can then be expanded in terms of the dressed 
states according to [$) — X]n=i *«('') I Xn(i"))j where 
the wavefunctions ^ — (^i, ^2, • • • , ^w+i)^ obeys the 



■f.9 T 



:^(P-A)2 
2to^ ' 



V 



* 



(28) 



with V being an external trapping potential. Here the 
vector potential A and scalar potential V are {N + 1) x 
{N + 1) matrices come from the spatial dependence of 
the atomic dressed states with elements 

A„,„ =: ifi,(x„(r)|Vxm(r)) , (29a) 

Vn,m = e«(r) Sn,m + {xn{r)\V (r)\xrn(r)) ■ (29b) 

If the off-diagonal elements of the matrices A and V 
are much smaller than the difference of the dressed 
atomic energies, we can safely neglect the off-diagonal 
contributions. In this regime, atoms in any one of the 
dressed states adiabatically evolve according to a sep- 
arate Hamiltonian with a U(l) Abelian gauge potential, 
like a charged particle moving in an electromagnetic field. 

The key point to generate non- Abelian gauge poten- 
tials is that some atomic dressed states, say the first 
q ones, form a degenerate (or nearly degenerate) man- 
ifold, leading to the failure of adiabatic approximation 
since some off-diagonal coupling can no longer be ig- 
nored. Suppose that these levels are well separated from 
the remaining ones, the full Hamiltonian can be pro- 
jected to this subspace for the reduced column vector 
^ ~ (^1, . . . , ^q) , which obeys the Schrodinger equa- 
tion ih-Sr^ ~ H^ with the Hamiltonian 



"dt 



Ira 



(30) 



where A and V being the truncated 9 x g matrices with 
elements defined in Eq. (29), and a scalar potential $ 
arises from the projection, with elements given by 






1 ^ 
-— V A„,,-A,^, 
2m ^-^ 



(31) 



l=q+l 



Here the effective U(q) non-Abelian gauge potential A 
also called the Mead-Berry connection is related to an 
effective magnetic field (or curvature) B: 



B = V X A + —A X A. 
in 



(32) 



Because the vector components of A do not necessarily 
commute, the term A x A does not vanish in general, 
even if the vector potential is uniform in space. In fact, 
this term reflects the non-Abelian character of the gauge 
potentials and appears only if g > 2. When q = 2 cor- 
responding to a tripod scheme discussed in the following 
section, under certain conditions the non-Abelian gauge 
structure is equivalent to a spin-orbit interaction. 
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Fig. 9 Schematic of atom-laser interaction, (a) Multi-pod con- 
figuration. An atomic state |e) is coupled to N different atomic 
states \gj) {j = 1, ..., A'') by N resonant laser fields. (b),(c) Atoms 
witii tripod-level configuration (N = 3) interacting witii three laser 
beams characterized by the Rabi frequencies f2i , f22 and ^3 . (d) , (e) 
Atoms with A-level configuration interacting with laser beams char- 
acterized by the Rabi frequencies Q'j, Qj a^cl a large detuning A. 

The above single particle physics of generating gauge 
fields is equivalent for bosonic and fermionic atoms. For 
a many-body system of bosons, such as a BEC, the U(l) 
Abelian gauge potential acts like a magnetic field and 
leads to a vortex lattice in the BEC [98,109]. The NIST 
group has recently realized the Abelian gauge poten- 
tial in a BEC of ^'^Rb atoms [110,111]. In their experi- 
ments, they used a gradient of the detuning of the laser 
frequency for the Raman transitions between different 
atomic ground sub-levels to generate an artificial gauge 
potential, somewhat different from the scheme we dis- 
cussed above. Besides, the NIST group has simulate 
an electric force acting on neutral atoms through the 
time dependence of an effective vector potential [112]. 
A many-body system of pseudo-spin- 1/2 bosons within 
the tripod scheme was first considered in Ref. [113], and 
the degenerate condensation ground state of the system 
was found in the presence of weak interactions. Further- 
more, in Refs. [114,115], a spin-orbit coupled spinor BEC 
was investigated and found to have two different phases 
in the ground state depending on the interatomic inter- 
actions, a single plane wave phase with a single dressed 
state and a standing wave with a nontrivial spin stripe. 
Most excitingly, the NIST group in a very recent exper- 
iment has reported the first realization of an effective 
U(2) non- Abelian gauge potential in neutral atoms - a 
spin-orbit coupled BEC [116]. 



B. Tunable Dirac-type equation for cold atoms 
with tripod- and A- level structure 

Specializing the above scheme of generating U(N) non- 
Ablelian gauge potential in neutral atoms, we show in 



this section two different systems, where a synthetic U(2) 
non- Abelian gauge field appears, can use to simulate the 
relativistic Dirac equation. 

Tripod-level configuration. — Let us first consider the 
adiabatic motion of atoms in x-y plane with each hav- 
ing a tripod-level structure in a laser field as shown 
in Fig. 9(b) and (c) [100,102]. The atoms in three 
lower levels |1), |2) and |3) are coupled with an ex- 
cited level |0) through three laser beams character- 
ized by the Rabi frequencies f^i = ilsinfle"'''^/-^^, 
^2 — ^sinOe^'^^ /V^ and fls = ficos^e"*"^, respectively, 
where fl = v^|rii|2 -t- [fial^ + [f^sl^ is the total Rabi fre- 
quency, and the mixing angle 6 defines the relative inten- 
sity. The Hamiltonian of a single atom reads 



H 



p2 

2m 



Hi, 



(33) 



where V = X)7"=i y?('')b)OI is the external trapping po- 
tential and the laser-atom interaction Hamiltonian Hint 
in the interaction representation is 

iJ„t = -?i(f]i|0)(l|-t-r!2|0)(2|+r!3|0)(3|)+H.c. (34) 

Diagonalizing the Hamiltonian Hint yields two degener- 
ate dark states of zero energy 



1 



|Di) = ^e-*"^ (e'""|l) - e-'""|2)), 
v2 



(35a) 



ID2) = ^e-*''^cos6'(e''^^|l)+e-"'^|2))-sin6l|3),(35b) 

as well as two bright state separated from the dark states 
by the energies zthft. If fl is sufficiently large compared 
to the two-photon detuning due to the laser mismatch 
and/or Doppler shift, then one can safely study only the 
internal state of an atom adiabatically evolves within the 
dark state manifold. In such a resonant coupling configu- 
ration, there is atomic dissipation due to the spontaneous 
emission of the excite state, thus the dark states acquire 
a finite lifetime. However, it can be up to a few seconds 
for cold atoms with typical velocity being of the order of 
a centimeter per second. The lifetime of the dark states 
can be further increased by using off-resonant coupling 
configuration as described below. 

In the dark state manifold, the atomic state vector 
1$) can be expanded in terms of the dark states ac- 
cording to 1$) = I]^^i *j(i")|£'j(r)), where *j(r) is 
the wave-function for the center-of- mass motion of the 
atom in the j-th dark state. Thus the dynamics of 
the atom is described by a two-component wavefunction 
^ = (\l/i, \['2)"^, which obeys the Schrddinger equation 
(30) with q = 2. The non- Abelian gauge potential A in 
the present configuration of the light field is of U(2) non- 
Abelian one. Substituting Eq. (35) into Eq. (29) and 
(31), one can work out the expression of potentials A, $ 
and^: 



A = Hk 



-e,.cos( 



—63; cos 9 
By cos^ 9 



(36a) 
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^^f n^K^ sin^ e/2 



m 

2,^2 „• 2 



ft^fi;^sin^(20)/8TO 



V 



Vi 







Fi cos2 6* + t/3 sin^ 



(36b) 
(36c) 



where the trapping potential is assumed to be the same 
for the first two atomic states, V\ = Vi. Choose Vx — Vj, = 
H'^k'^ sin 6 /2m, and thus the overah trapping potential 
simplifies to F + $ = V^il with I being the unit matrix. 
Furthermore, let the mixing angle 9 — Oq to be such that 
sin^ 9q = 2 cos 9o, giving cos ^o = \/2— 1. Thus the vector 
potential A takes a symmetric form 



A = Hk' { — ex<7x + Gj^Cr^) + hKQByl 



(37) 



in terms of the Pauli matrices a^ and a^, where k' = 
Kcos^o and kq = k(1 — cos^o)- Using a unitary 
transformation H' ^ W HU in Eq. (30) with U = 
exp(— iKoj/) exp (— i^cTj;), one obtains the Hamiltonian 
for the atomic motion 



2m 



(38) 



where a± ~ GxCTx + g^ctj, is the spin 1/2 operator in the 
xy plane. In the limit of low momenta, i.e., |P| <C ft/t', 
one can safely neglect the kinetic energy term in Eq. (38) 
and obtain an effective Dirac Hamiltonian 



H2D = ficok • (T_L + Vi + mcQ, 



(39) 



where cq = hn' /m is the effective light velocity for such 
quasi-relativistic particles. In this case, the atoms in the 
system behave as 2D massless relativistic particles. 

If we focus on the case with an extremely anisotropic 
potential, which is sufficiently strong along the trans- 
verse direction of the x axis, such that the original mas- 
sive atoms are actually confined in a one-dimensional 
guide along the x axis with the vector potential A^ = 
—hKC0s9ax- In this case, we get a ID Dirac equation 
described by [103,104] 



Hid = Ci,(JxPx + Izi^zi 



(40) 



up to an irrelevant constant, provided that the wave vec- 
tor of the atoms Px/ti <C ncosO. Comparing the original 
Dirac equation, here 7^ = h^K^ sin'* 9 /2m, is the effective 
rest energy and c^, = hn cos 9 /m is the effective light ve- 
locity, and thus the effective mass m^ = ^ tan^ 9 aiv? 9. 
Furthermore, if we take Vi — V3 = ft^K^ sin^ 9 /2m, as dis- 
cussed above, we are able to obtain a massless ID Dirac 
equation when 7^ in Eq. (40) vanishes. Note that the 
mass rnii, of the simulated Dirac particle is not the mass 
m of the cold atom itself and it is a remarkable feature 
that the mass term in the simulated Dirac-like equation 
can be controlled by the laser beams, so one can simulate 
a tunable Dirac equation with cold atoms. 

A-]evei conGguration. — We now turn to consider the 
adiabatic motion of atoms in y-z plane with each hav- 
ing a A-level structure interacting with laser beams as 
shown in Fig. 9(d) and (e) [105,117]. The ground 



states |1) and |2) are coupled to the excited state |3) 
through laser beams characterized respectively with the 
Rabi frequencies il[ — il' cos{Kyy)e~'^'^''^ and fl'2 = 
n' sm{Kyy)e'^''-'^-^\ where n' = ^\Q[\^ + \n'^\^. Such 
Rabi frequencies can be realized with two pairs of plane- 
wave lasers as shown in Fig. 9(e). For instance, il.'i can be 
realized with a pair of lasers with ^J7' exp[i{— k^zzL Kyy)], 
while il'2 is achieved with ^17' exp[i(±7r/2 — KzZ ± Kyy)], 
where Ky = k cos ip and n^ — n sin ip with n being the 
wave number of the laser beams and (p being the an- 
gle between the laser beams with the y axis. Such laser 
beams means that fl[ and Q2 are standing waves in the 
y direction but plane waves in the z direction, and both 
propagate opposite to the z direction but the phase of 
fl[ is TT ahead. The interaction Hamiltonian in Eq. (33) 
in this configuration is given by 






/iA|3)(3| 



2 
(^fi^.|3)(j| 



h.c), 



(41) 



where A is the detuning. Diagonalizing H[^^ 

yields a dark eigenstate \x^) = sin(«;j,2/)|l) + 
cos(Kj,y)|2) and two bright eigenstates \xi) = 
(sina) cos(Kj,y)|l) — (sina) sin(Kj,j/)|2) — cosQ;e^*'^^^|3), 
Ixf) = (cosa)cos(Kyy)|l) — (cosa) sin(Kj^?/)|2) + 
(sina)e~*'^''^|3), where a is determined by tana = 
(VA^ + Vt'"^ — A)/J7'. The corresponding eigenvalues are 
^^ == and EI2 = ft [A ± VA2 + A^''^\ /2, respectively. 
In the large detuning case |A| ^ i7', the two lower eigen- 
states Ix^) and Ixf ) span a near-degenerate subspace, 
and have negligible contribution from the state of E^ . 
Thus one can safely consider the adiabatic motion of 
atoms in the near-degenerate subspace, leading to a U(2) 
non-Abelian gauge potential with the vector and scalar 



potentials obtained as A 



/lK,£l' 1 



2A2 



(JyEz 



I blxiiu (J ■■(jCfj 



~2^' 



and $ = 



rv'n'^K^ 



4mA^ [cos(2i/?)cr^ + I] in this near- 
degenerate subspace with the basis {|x^), jxf )}■ Under 
this condition, we obtain an effective Hamiltonian 



Hh 



pI 



2m 



+ v'y(7yPy+v'z(7zPz+l'z<yz+VT, (42) 



where the parameters v'^ 



^z 

hn'- 



2mA2 ' 



and 



rz = T^^i^y - (1 + ^'V^^Wz] + ^+Vt. In the 
derivation, one has dropped an irrelevant constant and 
assumed that the trapping potentials Vi,2,3 = Vt are 
state-independent. The term related to azPz in Eq. (42) 
can be dropped approximately due to Vy ^ v'z- In fact, 
the dynamics of atoms in the z axis is still governed by the 
Schrodinger Hamiltonian since the kinetic energy term 
dominates, and thus the atoms can be well confined by 
an optical waveguide along y axis [103]. In the limit of 
low momenta, i.e., Py <C ft/tj,, one can neglect the kinetic 
energy term. To this step, one obtains the ID effective 
Hamiltonian 



H'lD = '"y'^yPy + iz^z + Vt, 



(43) 
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where v' is the effective speed of light in cold atoms and 
7^ is the effective rest energy, both of which are also 
controllable by the laser fields. 

In the end of this section, we note that the relativistic 
quasiparticles with total pseudospin 5 = 1, which have 
been described in Section II. A. 2, can also be simulated 
with the tetrapod (TV = 4 in Fig. 9(a)) scheme [118]. 
However, it looks like that the present multi-pod coupling 
scheme is not feasible to simulate the relativistic quasi- 
particles with higher (than 1) or even arbitrary pseu- 
dospin [118]. To overcome this limitation, one may de- 
sign other atom-light interaction schemes which remain 
further work, or one can adopt the proposal described 
in Ref. [119], where the massless Dirac fermions with ar- 
bitrary spin can be simulated by using multicomponent 
ultracold atoms in a two-dimensional square OL. 



C. Simulation of nonlinear Dirac equation 

In the previous section, we have shown that neutral cold 
atoms in a light-induced non-Abelian gauge filed can 
be used to simulate the 2D and ID relativistic Dirac 
equation with tunable parameters. Note that the above 
discussion is the single particle physics and is suitable 
for both noninteracting bosons and fermions. However, 
atomic interactions are natural and crucial in the atomic 
assemble, such as BECs. Besides, the artificial non- 
Ableian gauge potential has been realized in BECs [116]. 
Thus one may ask whether the relativistic quasiparticles 
can emergence in a many-body cold atom system. In this 
section, we shall address that in the presence of a laser- 
induced spin-orbit coupling an interacting ultracold BEC 
may acquire a quasi-relativistic character described by a 
nonlinear Dirac-like equation [105,120]. 

In the tripod scheme [120], one may choose another 
three laser beams to get a more symmetrical gauge poten- 
tial A ~ —hndyex- This may be achieved by employing 
three co-propagating lasers along the z axis, with con- 
stant [fia] and spatially dependent transversal profiles 
fii = [rial cos 0(a;,y)e*''^, Vt2 = [ris] sin^(a;, j/)e*''^ with 
a laser density modulation along x such that 0(r) — 
\i1kx. If the density modulation has however a polar 
symmetry on the xy plane, i.e., (/'(r) = \l2Kp (with 
p2 — 2;2 _|^ y2^^ then A = —hKayCp, which is the 2D 
case. Finally, setting state-dependent trapping potential 
Vi — y-i = fi^T — ?i^K^/2m with a detuning At, and 
V3 = — Vi — 2hAT, one obtains for both ID and 2D ar- 
rangements ^ + V — HAtctz- 

Interactions are characterized by s-wave scattering, 
and the scattering lengths of different internal states 
are in principle different, which may present collision- 
ally induced spin rotations. For simplicity of the dis- 
cussion, we consider that the interaction among the 
particles of different internal states are all identical, 
with an effective d-diniensional interacting strength g ~ 
4:TTh'^asN/{m\/2TTl'j_)^~'^, where a^ is scattering length, 
N is the particle number, and l± is the oscillator length 



associated to a harmonic vertical confinement. We also 
assume that the interaction energy is much smaller than 
hQ, such that we remain in degenerate subspace. Con- 
sidering a wave-packet with (p) = hko <^ Hk and mo- 
mentum width dp <C 2ftK, we can still safely neglect the 
p^ term. Within the Gross-Pitaevskii formalism, the in- 
teracting bosons in the degenerate subspace manifold are 
then effectively described by a NLDE as ihdt'^ = Hfqo^, 
where 



Hnd = — P • o- + HAtctz + .9*^ • *• 
m 



(44) 



Here the spin-orbit coupling term p-cr = Px<yy for the ID 
case, and for the 2D case p • cr = pa-y with p^ = Px + Py- 
For the A- level scheme [105], the requirements of 
weak atomic interactions and small momentum width 
for obtaining the NLDE are similar with those in tripod 
scheme, i.e., the interaction energy is much smaller than 
hfl' and dpy <C 2hKy. In this configuration, the dynamics 
of a spin-orbit coupled BEC are effectively described by 
the ID NLDE with the Hamihonian 



Hnd 



hv'aydy + 7>2 + g*^ • * + Vt- 



(45) 



Up to this, we have addressed that the dynamics of a 
BEC with a light-induced spin-orbit coupling can ef- 
fectively described by a ID or 2D NLDE under cer- 
tain conditions. Such quasi-relativistic BECs may be- 
come self-trapped and resemble the so-called chiral con- 
finement studied in the context of the massive Thirring 
model [120]. In Section IV. C, however, we will show that 
macroscopic relativistic tunneling effect can be detected 
in the A-scheme system. 

Before closing this section, we note that the NLDE may 
also be simulated with a BEC loading in the hexagonal 
optical lattice, proposed in Refs. [121-123]. It was sug- 
gested that the NLDE can be obtained by cooling bosons 
to form a condensation in the lowest Bloch band of a 
hexagonal optical lattice and then adiabatically trans- 
late the BEC to the Dirac point at the band edge by adi- 
abatically tuning the relative phases between the laser 
beams [121,122]. In this system, the interplay between 
nonlinear atomic interactions and Dirac dynamics may 
lead to rich localized excitations [122], including solitons, 
vortices, skyrmions, and half-quantum vortices. How- 
ever, it was recently argued that the tight-binding model 
to describe such an interacting boson system around the 
Dirac point is inadequate [123]. To support this point, it 
was showed that the atomic interaction of arbitrary small 
strength can completely deform the topological structure 
of the Bloch bands at the Dirac point [123]. 



IV. OBSERVATION OF SOME RELATIVISTIC 
EFFECTS 

Some relativistic effects related to the Dirac equation, 
have been predicted for a long time, the most well-known 
ones should be the Zitterbewcgung (ZB) effect and the 
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Klein tunneling (KT) . The concept of ZB was first intro- 
duced by Schrodinger [39], who calculated the resulting 
time dependence of the electron velocity and position in 
the framework of the Dirac equation. He found that, in 
addition to classical motion, a free relativistic electron ex- 
hibits a rapid periodic oscillations. The ZB effect is due 
to an interference between positive and negative energy 
state [38], and will disappears with time if an electron 
is represented by a wave packet since the interference 
has a finite lifetime. Though Schrodinger's idea stim- 
ulated numerous theoretical investigation, direct exper- 
imental observation seems impossible for free electrons 
since the predicted frequency fkoz — ^nieC^ — 1 MeV 
(here nie is the mass of an electron) and the amplitude 
Ac = h/mec ~ 3.86 x lO'^ A. 

Another exotic relativistic phenomena in the Dirac 
equation was pointed out by Klein, who first used the 
Dirac equation to study an electron scattering by a po 
tential step and found that there exists a nonzero trans 
mission probability even though the potential hight tend 
to infinity [40], in contrast to the scattering of a non 
relativistic particle. This unique scattering process ha 
attracted lots of interest over the past eighty years bu 
failed to be directly tested by elementary particles du 
to the requirements of currently unavailable electric fielc 
gradients, such as on the order of 2r7ieC^/eAc ^ 10^^ V/n 
for free electrons. 

In this section, we demonstrate that the famous ZI 
and KT can both be observed in cold atom systems 
where the dynamics of the quasiparticles are describe! 
justly by the relativistic Dirac equation. It is a remark 
able feature that the effective mass and the effective ligh 
speed are both controllable by the laser fields, whicl 
makes the observation region can be achieved in experi- 
ments. In principal, the OL systems discussed in Section 
II and the bulk atom systems without OLs discussed in 
Section III, are both able to exhibit the relativistic ef- 
fects, we however here focus on the latter for consistency 
but without loss of generality. We are also interested 
in the tunneling issues of a spin-orbit coupled BEC de- 
scribed by a NLDE, which may support the direct obser- 
vation of macroscopic KT under realistic conditions. 



A. ZitterbeweguDg 

The concept of ZB was first rooted in the Dirac equation, 
however, its present is not unique to the relativistic elec- 
trons, but rather than a generic feature of spinor systems 
with Dirac-type dispersion relationship. In solid state 
systems, the ZB can be viewed as a special kind of inter- 
band transitions with generation of virtual electron-hole 
pairs [124]. The energy gaps in these systems are usually 
of the order of meV and the oscillation frequency is typ- 
ically a few THz. Particularly, semiconductor quantum 
wires [125], graphene [126,127], superconductors [128], 
and photonic crystals [129,130], trapped ions [10,131] 
have been proposed as candidate systems for the obser- 



vation of ZB. Excitingly, a proof-of-principle quantum 
simulation of the ZB has been respectively performed in 
the recent experiments in a photonic crystal [131] and 
trapped ion system [10]. Since a Dirac-like Hamilto- 
nian can be created by cold atoms coupling with laser 
beams, atoms in tripod-scheme system undergo ZB as a 
result [104,132]. It is also possible to realize ZB with cold 
atoms in an Abelian vector potential [133], and even to 
control the amplitude and lifetime of ZB via an additional 
mirror oscillation [134]. The ZB can also be simulated 
with cold atoms in a driven [135] or titled OL [136]. In 
these cold atomic systems, typically, the frequency of ZB 
is on the order of a few MHz or kHz and the amplitude is 
comparable to the wavelength of the laser beams. Thus 
they are promising candidates for observing ZB. Here we 
give a complementary example to show the atomic ZB 
by using the A- scheme. 




Fig. 10 Zitterbewegung for wave packets of ^Li atomic gas with 
initial wave number fco = 0. (a) Ky = 5.0 X 10^ m^^, (b) Ky = 
1.0 X 10^ m^^. Other parameters arc m = 1.16 X 10^^®, l^/ft = 10 
kHz and (5; = 10 ^m, respectively. 

In the A-scheme, the dynamics of atoms are effectively 
described by the Dirac-like Hamiltonian (44). We con- 
sider the Gaussian wave packet as the initial wave func- 
tion 



*(y,0) = 



1 



VS. 



= e*'^"2'e' 



(y- 



iV^ 



C2 



(46) 



in the spatial space, where ]ci]^ (i = 1,2) normalized as 
kil + IC2I — 1 determine the initial population in states 
\Xi)i ^h ko and j/o a-re initial wave-packet width, initial 
averages of wave number and position, respectively. We 
suppose Si = -y/ft/TOWT, corresponding to the ground- 
state width of the trapping potential ^mcj'^y'^ in y-axis. 
In the momentum space, the initial wave packet is given 

by 



Hky.O) 



1 



*(y,0)e-''^««dy 



\/5kV^ 



C2 



(47) 



where 5k = 5j^ is the momentum spread. When t — Q, 
one turns off the trapping potential Vt, and after t time 
evolution governed by Dirac-type Hamiltonian H[jj with 
Vt = 0, the finial wave function is written as 
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*(y, t) = re-^K'^«P«+^>-)**(y, 0), 



(48) 



where T denotes the time ordering operator. It is 
straightforward to show its time evolution ^(y,i) — 
-^mky,t)e^>'yydky with 



<^{ky,t) 



_^-i{ky~ka)yo^~{ky-ka)^ /2Sl 



ci cos(a;fe<) " " ^^ ' sm(wfet) 

C2 cos(wfct) + " " ' ' sm(wfct) 



(49) 



where Wfc = \ {^U^Y + Wy^yY leads to population 
transfer between two spin states, and thus "^{y^t) is sen- 
sitive to the initial spinor components. If we choose 
ci — C2 — l/v2, the wave packet of atomic gas undergoes 
ZB [104,130]. To investigate atomic ZB, we calculate the 
expectation value of the center of mass, which is given 
by 



{y{t)) = i / dky(^\ky,t)dky^{ky,t) 



IK) 



hSkV 



dh 



sin'(^fcO (fc.,-fc„)= 



/si 



(50) 



The integral term stands for ZB which shows the oscil- 
lation of the center-of-mass motion. To have some ideas 
about ZB, we numerically calculate Eq. (50) and two 
examples are shown in Fig. 10. The reason for ZB is 
the separation of initial wave packet by spin and then 
interfere between the coupling components, which gives 
rise to oscillation of center-of-mass motion. ZB oscilla- 
tions for finite momentum spread damp out over time, 
and its lifetime can be increases by reducing momentum 
spread &k- In this model, the amplitude of ZB is pri- 
marily determined by the central momentum fco and the 
spin-orbit coupling strength v' . The smaller fco and the 
larger v'y, the larger amplitude of ZB. The amplitude of 
ZB here can reach a few /xm and it is around 7 /im when 
fco = 0, as shown in Fig. 10(b), which is larger than 
that in early proposals on observations of ZB with cold 
atoms [104,132-134]. Meanwhile, the atomic ZB can per- 
sist for several milliseconds with frequency being of the 
order of a few kHz. ZB of atomic gas with such large 
amplitude and long lifetime is detectable in current ex- 
periments. 



B. Klein tunneling 

The KT is originally referred to the step potential, how- 
ever, it has been extended to other kinds of potential 
barriers, leading to a general description that relativistic 
particles can penetrate through high and wide potential 
barriers without exponential damping expected in non- 
relativistic tunneling processes. Though this interesting 
relativistic scattering process can not be directly tested 
by elementary particles, some quasiparticles in certain 
systems may be described by effective relativistic wave 



equation and thus provide a platform to simulate the 
KT. Those including electrons in graphene [137], electro- 
magnetic waves in honeycomb photonic crystals [138], 
cold atoms in a driven [135] or bichromatic [139,140] 
optical lattice, confined stationary light [141], and trap 
ions [131], have been proposed to observe such rela- 
tivistic tunneling. Interestingly, in two very recent ex- 
periments, a proof-of-principle simulation of KT with 
trapped ions [11] and with ultracold atoms in a bichro- 
matic optical lattice [140] has been performed, respec- 
tively. We in the following review a proposal of directly 
observing the KT with cold atoms in a non-Abelian gauge 
field generated by using the A-scheme [105], as discussed 
in Section III.B, where we have obtained a Dirac-like 
equation (43) with tunable parameters. 
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Fig. 11 Schematic illustration of the system, (a) Atom with 
A-level configuration interacting with laser beams, (b) The con- 
figuration of the laser beams to realize a Dirac-like equation and 
an effective Gaussian (square)-shape potential induced by another 
laser beam. The atoms are confined in a ID waveguide along y axis 
and scattered by the potential. 

To have an intuitive physics picture, we first consider 
a single atom with energy E scattered by a square poten- 
tial with the width L and potential hight Vg, as shown 
in Fig. 11(b). Such an effective square potential can 
be experimentally formed by a laser beam with flat-top 
profile [142]. The transmission coefficient Td for the so- 
called KT regime is given by 



Td= [l + {ii-Ti-^fsin\aL)/A] 



(51) 



where rj 



{V.-E+^'^)(E+^'^) ^ 



(^-£-7^) 



^Jh. 



(V,-E+i'^)l 

Compared with the familiar property in the textbook of 
non-relativistic quantum mechanics that the transmis- 
sion coefficient is a mono exponential decreasing as a 
function of the width L or the potential height Vs , a dis- 
tinguished different feature within this relativistic tun- 
neling region is that the tunneling amplitude can be an 
oscillation function of Vs or L even when the kinetic en- 
ergy of the incident particle is less than the height of the 
square barrier potential. This relativistic effect can be at- 
tributed to the fact that the incident particle in positive 
energy state can propagate inside the barrier by occupy- 
ing a negative energy state, which is also a plane wave 
aligned in energy with that of the particle continuum 
outside. Matching between positive and negative energy 
sates across the barrier leads to the high-probability tun- 
neling. Furthermore, the resonant transmission occurs 
for perfect matching where the potential width equates 
integral multiple of half-wavelength of the negative en- 
ergy state, corresponding to aL = nir (n = 1, 2, • • •) in 
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the formula of Tp. We take the atoms of ^Li with mass 
m = 1.16 X 10~^^ kg as an example. If we choose the fol- 
lowing typical experimental parameters: the wave num- 
bers ka = 1.1 X lO*' m-\ Ky = 10^ m-\ k^ = 0.8 x lO"^ 
m~^, the Rabi frequency ft = 10^ Hz and the detuning 
A = 10^ Hz, we can find that the Klein regime cor- 
responds to the required potential height Vs > 0.162 
?ixMHz, which can be easily achieved in experiments. 
So we have demonstrated from a simple example that it 
is feasible to observe the KT with cold atoms. 

As for a practical experiment we are required to re- 
lease two conditions: the trajectory of a single atom is 
hard to detect, and it is much easier in experiments to 
measure the density evolution of an ensemble of nonintcr- 
acting atoms. Compared with a perfect square potential 
barriers achieved with a flat-top laser beams, a Gaus- 
sian potential barrier V^{y) — Vce"^ Z"' , where Vg is 
the height and cr characterizes the corresponding spa- 
tial variance, is much easier to be generated. However, 
the conditions of resonant transmission varies with the 
velocity and the width of the potential, and thus both 
the ensemble of atoms and the Gaussian potential may 
smooth the oscillations in the transmission coefficient. 

So we now turn to address an ensemble of atoms that 
are scattered by the Gaussian potential. We assume that 
the ensemble of atoms are initially trapped in a harmonic 
trap which moves along the y axis with the wave number 
ka- At the beginning t — 0, the center of the harmonic 
trap locates a.t y = —d, and the center of the Gaussian 
potential barrier is at y = 0, as shown in Fig. 11. In this 
case the number density and the number wave number 
of the atomic ensemble are characterized by Gaussian 
distribution. The trap is turned off at i = and then we 
calculate the evolution of the density profile of the atomic 
gas after a long enough time for scattering. We define 
the average transmission coefficient for an ensemble of 
Na noninteracting atoms as 



(T) 






Na 



(52) 



where T{ki) denotes the transmission coefficient of the 
atom i scattered by the Gaussian potential, as a function 
of random wave number ki described by the Gaussian 
distribution ki ~ N{ka,cff.). The analytical expression 
of r(fci) is absent, however, we here show an efficient 
method to numerically solve it based on the transfer ma- 
trix methods. The numerical procedures are outlined as 
follow. One first cuts the Gaussian potential into spa- 
tially finite range y € [— ^cyc], where the cutoff posi- 
tion j/c should be chosen to guarantee that the potential 
height outside the range is low enough to be transpar- 
ent for the atoms, i.e., V^{yc) ^ E.Vq. Secondly, one 
equally divides this range into n spindly segments and 
each segment can be considered as a square potential if 
n is large enough. Thus the Gaussian potential can be 
approximately viewed as a sequence of connective small 
square potential barriers, and the transmission coefficient 



T(ki) « l/|mii|^, where mn is the first element of the 
whole transfer matrix M = M„M„_i • • • M,- • • • M2M1. 
Here Mj denotes the transfer matrix of the j-th square 
potential barrier, whose explicit elements can be found in 
Ref. [103]. Note that this numerical calculation scheme 
recovers the non-relativistic scattering governed by the 
Schrodinger equation. 
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Fig. 12 Klein tunneling of an ensemble of Na = 10* noninteract- 
ing atoms, (a) (T) as a function of potential width a with f2p = 0.8 
MHz. (b) (T) as a function of potential hight Vg with cr = 10 fiva. 
The black solid line, blue dashed line and red dotted line in both 
(a) and (b) correspond to the cases of o-f^ = 0.5, 1.0, 2.0ka, respec- 
tively. Other parameters are ka = 1.1 x lO'' m~^, tzy = 10^ m~^, 
K^ = 0.8 X 10^ m-i, frequency U = 10^ Hz, and A = lO" Hz. 

Figure 12 shows the average transmission coefficient 
(T) of an ensemble of iVa = 10'* noninteracting atoms as a 
function of the width a and the height Vg of the Gaussian 
potential. Similar with the results for the single atom, 
the tunneling oscillation still survives after the average 
of the random distribution of the velocities for the Gaus- 
sian potential. The amplitude of oscillation decreases 
with the increasing of the wave number variance. Such 
KT features can be attributed to the fact that the tun- 
neling of an ensemble of noninteracting atoms is equiva- 
lent to the superposition of every single-atom tunneling 
with different oscillation period and amplitude, leading 
to smoothing the oscillation property. The large veloc- 
ity variance, the more obvious smoothing effects. The 
oscillation period of single-atom tunneling is insensitive 
to the potential hight in contract to the potential width, 
and thus almost no decay of oscillation in Fig. 12(b). 
In particular the interesting KT can be observed under 
realistic conditions. For instance, the difference between 
the nearest peaks of (T) in spatial dimension are in the 
region 3.0 — 6.0 //m, and the frequency is about 0.2 MHz. 
Both of them in the oscillation are detectable within the 
current technology. 



C. Macroscopic Klein tunneling 

Although it is possible to observe the KT at a relatively 
high temperature, the phenomena is clearer in low tem- 
perature. In low temperature the bosonic atoms may 
form the BEC. Moreover, the gauge field has recently 
been generated in a BEC by the NIST group [111,116], 
thus it deserves to study whether the KT is still de- 
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tectable with a BEC. Surprisingly we will illustrate below 
that the KT of a spin-orbit coupled BEC may be observed 
very clearly [105]. 

The single-atom dispersion is characterized by two 
branches E±{ky) = ±(7^^ -I- h^v'^kyY/'^, where the lower 
(up) branch represents the negative (positive) energy 
state. One can prepare an initial BEC with a designated 
mode fco at the positive or negative energy branch. The 
two branches allow us to investigate a more fruitful tun- 
neling problem: there are four classes of the scattering 
which describe the wave function ^^ (m(= ±)) scattered 
by the potential vVb^ with v donating a barrier [v = +) 
or a potential well {v = — ), as shown in Fig. 13(a). We 
assume that the BEC is initially trapped in a harmonic 
trap which moves along the y axis, thus we may choose 
the initial wave function of the BEC as 



*m(2/:0) 



1 



vW^ 



gVfeo!;g-(a-(-d)^/2/^0 



(53) 



where Iq is the width, fco is the central wave num- 
ber of the wave-packet and the spinor cj)^ are defined 
as (j)+ — (zcos^, — sin^)^, (f>- — (— isin^,cos^)^ with 
^ — h arctan(/iWj,fco/7^) and T as the transposition of ma- 
trix. Equation (53) describes the Gaussian wave packet 
with the central velocity h{Ky + iJ,ko)/m moving along 
y-axis. After the evolution governed by Dirac-type (45) 
with time t and replacing Vr by i'V^{y), the finial wave 
function becomes 



* 



i^iy,t) =Texp ( -- 



H'jyjjdt 



*m(2/-0)- 



(54) 



We have numerically calculated ^/^(y, t) and as an exam- 
ple being shown in Fig. 13(b), the stationary solution of 
the tunneling is always found within a few milliseconds, 
which is a small time scale comparing with the lifetime of 
BECs. After tunneling, the incident wave packet divides 
into the left- and right-traveling wave packets and only 
the latter one is on the transmission side of the barrier. 
Thus we can define the transmission coefficient of the 
incident wave packet ^^(y,0) scattering by a potential 
lyVb'^ as 



T — 



'^i{y,T)^t.{y,T)dy, 



(55) 



where r (being slightly larger than cI/vq) represents 
a time that the reflected and transmitted wave pack- 
ets are sufficient away from the Gaussian potential. 
One can directly measure the transmission coefficient 
in Eq. (55) since the spatial density distribution 
Pfj.{y,T) = |\l/^(y,r)p can be detected using absorption 
imaging [143]. We first look into the tunneling phenom- 
ena for a non-interacting BEC {g = in Eq. (45)) real- 
ized by Feshbach resonance [143], and then briefly discuss 
the effects of the interaction between the atoms. 

We plot the transmission coefficient r_|_+ as a function 
of the potential height Vq and width a in Fig. 14(a,b) 
with the practical experimental parameters. It is inter- 



esting to note that the transmission coefficient decreases 
exponentially to zero with Vg when Vq < Vq , while if 
we further increase the potential height to the Klein re- 
gion Vg > Vq , the transmission coefficient will increase 
and then be an oscillating function in the Klein region 
Vq > Vq , being similar to the results of Eq. (51) for the 
which describes the transmission of the square potential 
barrier. Here the critical value of the potential height 
may approximately be estimated using the square poten- 
tial barrier with V^ = ^(fco) + 7z « 0.09 MHz. Besides, 

there are two identities T^^ = T and T ^ = r_| since 

Eq. (45) with g = is invariant under the charge con- 
jugation [144], the former as an example is confirmed in 
Fig. 14(a). Although the amplitude of tunneling oscilla- 
tion is less than the unit as comparing with the tunneling 
of single atom, the amplitude of tunneling oscillation can 
be more than 0.5 and meanwhile the period can be a few 
micrometers as shown in Fig. 14(b), which is experimen- 
tally detectable. 

The most exotic feature induced by the relativistic ef- 
fects is that, the BEC with negative energy can almost 
completely transmit a wide Gaussian potential barrier, as 
shown in Fig. 14(c). This phenomena can be attributed 
to the fact that the incident BEC in a negative energy 
level can be considered as a macroscopic 'anti-BEC Al- 
ternatively, it can be understood that such scattering fea- 
ture is actually equivalent to that of a BEC of positive 
energy scattered by a Gaussian potential well because of 
T |_ = T_| We also calculate the transmission coeffi- 
cient for the central mode of the wave packet, as shown 
in the insert of Fig. 14(c), which further confirms that a 
wide enough Gaussian potential well is transparent. The 
reason lies in the fact that, in contract to square potential 
wells, the Gaussian potential wells are smooth (without 
any energy jump) in the whole space, and may support 
adiabatic motions of wave packets in the large width lim- 
itation. 




Fig. 13 Schematic representation of the scattering of a BEC. 

(a) A schematic diagram shows four kinds of scattering events. 

(b) Normalized density distribution in a scattering process at time 
t = 0, 1.5 and 2.0 ms. The peaks at j; = are the Gaussian barriers. 

The exotic tunneling property exhibited in Fig. 14(c) 
is an intrinsic relativistic and macroscopic quantum phe- 
nomenon that can not be explained with an incoherent 
ensemble average of a large number of atoms. To clarify 
this point we also calculate numerically the average trans- 
mission coefficient for an ensemble of 10^ noninteracting 
atoms as (T) . Here we choose the wave number distribu- 
tion to be the same Gaussian distribution as that in Eq. 
(53), i.e., ki ~ A^(fco,(T^) with the variance tJfc = 1/Iq. 
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The average coefRcient (T) is shown in the inset of Fig. 
14(c), which is almost the same as that of a single atom 

since Cfc is small. The differences between (T) and T ^ in 

Fig. 14(c) clearly demonstrate that the tunneling of BEC 
is not equivalent to an ensemble average of the individual 
atoms with the same distribution of wave number. 



Before ending this section, we make two additional 
comments, (i) The quadratic term of the momentum 
in Eq. (42) has been neglected in the derivation of the 
Dirac Eq. (45). To judge the feasibility of this approx- 
imation, the transmission coefficients T-f + with or with- 
out the quadratic term are compared in Fig. 14(b). It is 
shown that the quadratic term leads to merely a slight 
left-shift of the tunneling peaks. This phenomenon can 
be interpreted by the fact that the wavelength of the 
'anti-particle' consisting of the BEC inside the barrier 
decreases slightly in the presence of the additional low 
kinetic energy. This result verifies that the approxima- 
tion which leads to the effective Dirac equation is well 
satisfied, (ii) We have also calculated the transmission 
coefficient for the weak atomic interaction in BECs in 
Fig. 14(b), which shows that the effect of the weak inter- 
action is little and smooths merely the tunneling oscilla- 
tion slightly. Therefore the exotic tunneling phenomena 
addressed here survive in the case of weak interaction 
between atoms. 
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Fig. 14 Klein tunneling of BECs. (a) T-|-+ as a function of the 

potential height Vq with cr = 5 ^m. The relation Tf + = T 

is confirmed, (b) T-|.+ as a function of potential width cr with 
Vq = 0.2 MHz. The tunnelings of a BEC with the classic kinetic 
energy term and the weak atomic interaction {N = 2x 10'', l± = 1.4 
^m and as = 5ao with ag being the Bohr radius) are also depicted. 

(c) T ^ as a function of potential width a with Vq = 0.2 MHz. 

In the insert, the transmission coefficient of a plane wave at bf the 
central mode fco is shown as the line with labels of triangle, while 
the average transmission coefficient (T) of 10* atoms is shown as 
the solid line. The other parameters: Iq = 10 /im, fco = 5.5 X 10^ 
m-l, 7^/^! = 30 kHz, and d = 4(/o + o"). 



V. CONCLUSIONS AND PERSPECTIVES 

In summary, we have reviewed recent progress on quan- 
tum simulation of the relativistic Dirac equation by us- 
ing ultracold neutral atoms. As we have shown, ultra- 
cold atoms in a wide range of systems, including the 
optical lattices with designated structures and different 
light-induced gauge fields, can behave as the relativistic 
particles under certain conditions. Comparing with other 
relativistic systems such as the graphene, cold atom sys- 
tems offer rather more degrees of freedom to the relativis- 
tic quasiparticles, which can exhibit from one-dimension 
to three-dimension, and even from massive to massless. 
In other words, one can simulate a highly tunable Dirac 
equation with cold atoms by well dressing laser beams. 
Furthermore, the controllable atomic interactions and 
disorders in these systems provide a unique platform to 
investigate the relativistic dynamics and scattering pro- 
cesses such as the celebrated Zitterbewegung effect and 
Klein tunneling as well as its extended macroscopic ver- 
sion. 

Remarkably, some crucial experimental techniques for 
cold atoms used to be Dirac equation simulators have 
been recently demonstrated. The important progress in- 
cludes addressing atomic gases in a hexagonal optical lat- 
tice [69], and generating an effective U(2) non-Abelian 
gauge potential in a BEC [116]. The latter technique is 
principally available for fermionic atoms, which may pro- 
vide an opportunity for the realization of the long-sought 
Majorana fermions (see Refs. [145,146] and the referred 
works therein) in ultracold atomic superfluid [147-149]. 
The Majorana fermions described by the Majorana equa- 
tion are distinguished from the Dirac fermions by the 
most exotic property that Majorana fermions are their 
own antiparticles. The Majorana fermions with the fea- 
ture of non-Abelian statistics has raised significant inter- 
est as candidates for the realization of topological quan- 
tum computation [150]. To achieve this goal, we will still 
face the arduous task of manipulating Majorana fermions 
in cold atomic superfluid as that in other proposed sys- 
tems [151]. Another important issue beyond the hunt 
for Majonana fermions should be the quantum simula- 
tion of Majoranna equation [152], which is central to the 
recent work not only in particle physics, supersymmetry 
and dark matter, but also some exotic states of ordinary 
matter. If succeed, we are allowed to explore some new 
exotic quantum relativistic processes that may go beyond 
ordinary Dirac quantum mechanics. 

In the present proposals of quantum simulation of the 
Dirac equation, the atomic interactions are less consid- 
ered due to the complexity such as the nonlinear inter- 
acting for bosonic atoms. Thus, it is of interest to ask 
whether and how the whole physics picture of relativistic 
dynamics in the present of controllable nonlinearity, that 
is an open question in the original Dirac theory, can be 
learned by cold atom simulators [105,120]. Besides, com- 
bining with the controllable laser-formed disorder poten- 
tial techniques, we may also be able to investigate the 
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Anderson localization of relativistic particles [103], es- 
pecially the relativistic version of competition between 
the interaction and disorder, which is still a controversial 
problem even for the non-relativistic particles. Further- 
more, shall we be able to realize the quantum simulation 
of interacting relativistic quantum field theories [153] or 
find some other exotic high-energy phenomena in these 
oppositely lowest-temperature systems among the uni- 
verse? 
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